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ABSTRACT 

We present a novel implementation of Smoothed Particle Hydrodynamics (SPHS) 
that uses the spatial derivative of the velocity divergence as a higher order dissipation 
switch. Our switch - which is second order accurate - detects flow convergence before it 
occurs. If particle trajectories are going to cross, we switch on the usual SPH artificial 
viscosity, as well as conservative dissipation in all advected fluid quantities (for exam- 
ple, the entropy). The viscosity and dissipation terms (that are numerical errors) are 
designed to ensure that all fluid quantities remain single-valued as particles approach 
one another, to respect conservation laws, and to vanish on a given physical scale as 
the resolution is increased. SPHS alleviates a number of known problems with 'clas- 
sic' SPH, successfully resolving mixing, and recovering numerical convergence with 
increasing resolution. An additional key advantage is that - treating the particle mass 
similarly to the entropy - we are able to use multimass particles, giving significantly 
improved control over the refinement strategy. We present a wide range of code tests 
including the Sod shock tube, Sedov-Taylor blast wave, Kelvin-Hclmholtz Instability, 
the 'blob test', and some convergence tests. Our method performs well on all tests, 
giving good agreement with analytic expectations. 

Key words: Multiphase Smoothed Particle Hydrodynamics, Numerical methods, 
Monte-Carlo methods. 



1 INTRODUCTION 

Smoothed Particle Hydrodynamics (SPH) is now widely 
used in almost all areas of theoretical astrophysics dGin- 



gold fc Monaghan 1977 Lucy 1977 Monaghan 19921). Its 



popularity has been largely driven by its Lagrangian nature 
that makes it manifestly Galilean invariant and geometry- 
free; its ease of implementation; and the fact that it couples 
naturally to tree-gravity solvers that are currently the most 
efficient method for solving gravity ( Monaghan||1992 



Price 



2005) |Rosswog|2009| [Springel|2010b| |Dehnen||2000| |Dehnen 
fc Read||2011| ). 

There are many different flavours of SPH used in the 
literature reflecting the above broad range of applications. 
The most common - that we shall call 'classic' SPH - is the 
fully conservative SPH implemented in the standard release 
of the GADGET-2 code ( |Springel fc Hernquist]|2002| |Springel 
2005f 1 Although classic SPH remains a powerful numeri- 



cal tool for solving the fluid equations, it suffers from slow 



et al.|2010||Springel|2010a| 



numerical convergence (Springcl 2010b), and a spurious sur- 



face tension at phase boundaries that inhibits fluid mixing 
(see e.g. Morris||1996a| |Dilts 



1999 



Agertz et al.||2007| |Wadsley et ai.||2008HPrice||2008| |Reac' 



Ritchie & Thomas 2001 



In recent work, we demonstrated that mixing in classic 



SPH fails for two distinct reasons (|Read et al.|2010| RHA10) 



The first is a leading order error in the momentum equation, 
previously identified by |Morris| (|1996b[) and |Dilts | fl!999|), 



that we called |Eo|. This can grow by orders of magnitude 



at flow boundaries, delaying the onset of in stabilitie s. The 
second is a pressure discontin uity at flow boundaries, previ- 
ously identified by|Ritchie fc Thomasj (|2001[), |Price| (|2008| 
j2008l~ 



and Wadslcy et al. (2008), that we called the local mixing 
instability (LMl|^] This leads to a large force error which 
manifests as a spurious surface tension. Both problems must 



* E-mail: justin.inglis.read@gmail.com 

1 Slightly different implementations of this algorithm are also 
used in the literature, for example in the Gasoline code | Wadsley | 



|et al.|2004| and the Hydra code Couchman et al. 1995] These are 

sufficiently similar to also be called 'classic' SPH. 

2 It is an instability since, even if we start in pressure equilibrium, 
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be solved in order for mixing between fluids of different den- 
sity or entropy to proceed correctly (see also fj3] and s|4] in 
this paper). 

In RHA10 we presented some simple proof-of-concept 
solutions to both of these problems. We cured the LMI by 
using a weighted density estimate first proposed by |Ritch"ie| 
& Thomas (2001), and we showed that |Eo| can be made ar- 



bitrarily small by brute-force so long as the method is stable 
to large neighbour number (this required introducing some 
new kernels). However, our resulting Optimised Smoothed 
Particle Hydrodynamics (OSPH) method required a neigh- 
bour number that scales linearly with the density contrast 
on the kernel scale. The OSPH pressure estimator is also bi- 
ased in regions of the flow where entropy gradients are large. 
This leads to poor performance in strong blast wave tests 
(we demonstrate this in Appendix |B|) . 

The above problems with SPH have led to a wel- 
come proliferation of new Lagrangian or pseudo-Lagrangian 
techniques in the literature, including a moving-mesh code 
(Springcl 2010a), flux-based particle metho ds (|Gaburov fc | 
Nitadori|2010f7S PH using a Riemann solver (|Inutsuka|2002[ 



Cha et al. 2010 



20111, and SPH using a 



Murante et al. 
Voronoi tessellation for the densities ( |HeB fc Spri ngcl 2010) 
It has also led to an exploration of improved flavours of 
SPH that add additional dissipation terms to mitigate the 
surface tension effect, and use switches to reduce the dissi- 
pation away from flow boundaries and s hock^j (e.g. Price 



2008} IKawata et al.|[2009| |Rosswog||20lol |Cullen fc Dehnen 



2010| |Price|20TT i" 

In this paper, we present a new flavour of SPH - SPHS 
- that has the mixing performance of OSPH, but does not 
introduce prohibitive numerical cost. As in OSPH, we use 
a larger than normal neighbour number with a correspond- 
ingly higher order and stable kernel to reduce the force er- 
rors. However, instead of the expensive OSPH pressure es- 
timator, we introduce a higher order dissipation switch to 
ensure that all fluid quantities are smooth by construction. 
We show that these simple changes to the SPH algorithm 
lead to converged results with increasing resolution, and ex- 
cellent performance across a wide range of code tests. Our 
dissipation switch also allows us to use multimass SPH parti- 
cles. SPHS is useful for any astrophysics application involv- 
ing multiphase fluid flow (e.g. resolving the ISM in galaxy 
discs), or where the use of multimass particles would be ad- 
vantageous. 

This paper is organised as follows. In Sj2] S|3]and [j4]we 
present the SPHS method. In S|5] we discuss our timestep cri- 
teria and multi-stepping scheme. In S|6| we describe our im- 
plementation of SPHS in the GADGET-2 code (|Springel|2005 1 



In SjT] we present a suite of tests for our new method that 
demonstrate that it can successfully model shocks, bound- 
ary instabilities, and shear flows. We also check that it con- 
serves momentum, energy and mass and discuss the numer- 
ical performance of the code. Finally, in S|8] we present our 
conclusions. 



an infinitesimal perturbation will cause a pressure discontinuity 
to form. 

3 Actually, some of these SPH flavours have been in use in the 



2 THE SPHS EQUATIONS OF MOTION 

In this paper, we consider solving the Euler equations in the 
absence of sinks or sources in the Lagrangian 'entropy form' 
( |Springel fc Hernquist|[2fJ02| : 

dp „ 
dt =^ V - V 

dv _ _ VP 
dt p 

A — const. (3) 



(1) 

(2) 



closed by the ideal gas equation of state: 
P = A{s)p y 



(4) 



where 7, p, v and A are the adiabatic index, density, velocity 
and specific 'entropy function' of the flow, respectively. The 
function A(s) is a monotonic function of the specific entropy 
s. For adiabatic flow in the absence of sinks or sources, A 
is conserved. Thus equation [3] implicitly solves the energy 
equation. If required, the specific internal energy can be cal- 
culated from A and p as: 



Ap 



,7-1 



7- 



(5) 



Note that often A is referred to as the 'entropy' when really 
it is a monotonic function of the specific entropy. From here 
on we will adopt this convention also. 

We use the discrete form of the above equations as in 
RHA10: 



E 

3 



m-jWij 



dt 

Pi = A iP l 



N 

E 

3 



PrP] 



[Pi + pa v t w l3 



(6) 



(7) 



(8) 



where m« is the mass of particle i; r»j = Vj — r^; Wij 
. [Wij (hi) + Wij (hj )] ; and W is a is a symmetric kernel that 
obeys the normalisation condition: 



WQr- r'|,/i)dV 



and the property (for smoothing length h): 



lim W(\r- 

h^0 



r'|, h) = S(\r — r' 



(9) 



(10) 



Note that we do not explicitly solve the continuity equation 
nor the energy equation. The continuity equation is implic- 
itly solved by equation [6] since its time derivative satisfies 
a discrete form of equation [l] (see e.g. Price|[2005 1 . The en- 
ergy equation is implicitly solved by advecting the entropy 
function Ai = const, along with the particles (Springel & 
Hernquist|2002"l > 



We use a variable smoothing length hi as in |Springet] 
& Hernquist (20021 that is adjusted to obey the following 



constraint equation: 



4tt 



h, m = N n 



with 



literature for quite some time (see e.g. Morris & Monaghan 1997 1 



N 

= E^ 



(ii) 
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where N n is the typical neighbour number (the number of 
particles inside the smoothing kernel, W). The above con- 
straint equation gives fixed mass inside the kernel if particle 
masses are all equal. 

The above equations of motion manifestly conserve mo- 
mentum, mass and entropy. They do not manifestly conserve 
energy, but the energy conservation is still extremely good 
as we will show in §7.7| A fully conservative form of SPH 
can be constructed by r eplacing equation [7| with eq uation 
Al (see Appendix |A| and S pringel fc Hernquist|2002 1. How- 
ever, as shown in RHA10, this leads to a larger truncation 
error in the momentum equation. In Appendix [Aj we show 
that - for the test problems presented in this paper - the 
fully conservative form gives only a modest improvement in 
energy conservation while introducing significantly more dif- 
fusion for multiphase test problems. For this reason, we use 
the above set of equations as our default choice for SPHS. 

So far, the above equations are very similar to classic 
SPH and thus will suffer from both the |Eo| error and the 
LMI problems described in !jl We now address each of these 
problems in turn in sections ! 3] and Q 
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3 ERRORS & CONVERGENCE 

The first problem with classic SPH is the |Eo| error in the 
momentum equation (RHA10). While this is minimised by 
using equation [7J the error is still present in SPHS. To see 
this, we can Taylor expand Pj about Pi in equation [7] to 
obtairjf] 



hi p.. 



E- 

3 



Pj 



+ 0(h) 



(12) 



where Vi is a matrix that approximates the identity matrix, 
and Vf = hi\7 is a dimensionless gradient operator. The left 
term in equation [12] defines the dimensionless En , error: 



En 



N 
3 



Pj 



(13) 



Taking the limit of infinite kernel sampling (and equating 
rrij/pj with a volume element dV), we see that En,j is a 
discrete approximation to the volume integral: 



En 



dVV x W = 



(14) 



which is zero because V^W is antisymmetric 

Although Eo.i should be approxi mate ly zero, it is prob- 
lematic because it appears in equation 



12 



mal convergence then requires that En,j shrinks faster than 
hi. This can be tricky to ensure and depends intimately 
upon the choice of kernel W employed. A popular choice is 
the cubic spline (CS) kernel: 



W 



irh 3 




< x sC \ 
\ < x SC 1 
otherwise 



(15) 



Figure 1. The CS (black), CT (red) and HOCT4 (blue) kernels 
and their first derivatives (dotted lines). The vertical lines mark 
the half mass radii for each kernel. For the CS and CT kernels 
the half mass radii overlap on this plot. 



where x = r/h and, as written above, h defines the kernel 
edge not its resolving power. (The two are not the same as 
can be readily understood by considering a Gaussian kernel. 
This has an infinite edge, but a resolving power given by ~ 
its scale length.) 

Now, it is tempting to increase the kernel sampling sim- 
ply by stretching h for the CS kernel. However, this is a bad 
idea for two reasons. Firstly, it introduces bias into the den- 
sity estimate, spoiling convergence. Secondly, the CS kernel 
is not stable to large neighbour number. As h is increased, 
the particles clump on the kernel scale and the sampling is 
not significantly improved. For these reasons, in RHA10 we 
proposed a new class of kernels that can be used to achieve 
convergence. The lowest order of these was the CT kernel: 



-12a + 18a 2 ) x + /3 <x^a 



W 



N 
J? 



1 - 
2(1 




6a; 2 + 6x 3 



a < x < j 
| < x < 1 
otherwise 



(16) 



at order h i \ For- where /3 = 1 + 6a 2 - 12a 3 , N = 8/[ir (6.4a 5 - 16a 6 + l)]; 



and a = 1/3. This has spatial resolution similar to the CS 
kernel but is stable to larger neighbour numbers. The next 
highest order was the HOCT4 kernel: 



W = 



N 
J? 



Px + Q 

(1 - x) A + (a-x) 
(1 - x) 4 + {a-x) 
{1-xf 




< X ^ K 
K < X < P 

P < x < a 

X < 1 

otherwise 



(17) 



4 Note that this assumes that the pressures are smooth and there- 
fore differentiable. In classic SPH, this is not guaranteed. How- 
ever, in SPHS, we add dissipation terms to ensure that this is the 
case. We discuss these in detail in ij4] 



with N = 6.515, P = -2.15, Q = 0.981, a = 0.75, = 0.5 
and k = 0.214. The CS, CT and HOCT4 kernels and their 
first derivatives are shown in Figure [l] 
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Table 1. Density errors for a selection of SPH kernels applied 
to a constant density box. The columns give: kernel type; neigh- 
bour number; lattice configuration (glass or simple cubic); and 
the median/5%/95% recovered density to two significant figures 
(the true density is ptrue = 1-00). 



Kernel 


N n 


Lattice 


P (5%) 


p (median) 


p (95%) 


CS 


32 


simple 


1.00 


1.00 


1.00 


CS 


128 


simple 


1.00 


1.00 


1.00 


CS 


32 


glass 


1.02 


1.01 


1.01 


CS 


128 


glass 


1.01 


1.00 


0.996 


CT 


32 


simple 


1.07 


1.07 


1.07 


CT 


128 


simple 


1.02 


1.02 


1.02 


CT 


32 


glass 


1.08 


1.07 


1.06 


CT 


128 


glass 


1.02 


1.02 


1.01 


HOCT4 


442 


simple 


1.01 


1.01 


1.01 


HOCT4 


442 


glass 


1.02 


1.01 


1.00 



We demonstrated in RHA10 that the HOCT4 kernel is 
stable for 442 neighbours on a lattice, while having a similar 
spatial resolution to the CT or CS kernel with 128 neigh- 
bours. This can be partially understood just from the half 
mass radii of these two kernels (Figure [TJ. If we assume that 
the resolution scale is the half mass radius, then 442 neigh- 
bours for the HOCT4 kernel is equivalent to ~ 240 neigh- 
bours for the CS kernel. However, in RHA10, we asses the 
spatial resolution more carefully by comparing instead the 
ability for these two kernels to correctly resolve the sound 
speed of linear waves. This is what leads to the conclusion 
that the HOCT4 kernel with 442 neighbours has similar re- 
solving power to the CS with 128. Compared to 'classic' SPH 
with 42 neighbours, the HOCT4 kernel with 442 neighbours 
has, therefore, a poorer spatial resolution of only a factor 
~ (128/42) 1/3 ~ 1.5. We will use the HOCT4 kernel with 
442 neighbours for our default SPHS scheme. The CT ker- 
nel with 128 neighbours will be used only for convergence 
testing. 



As pointed out by Price (20101, the CT and HOCT4 



kernels have slightly larger density error than the more stan- 
dard CS kernel. However, the effect is small (see Table [T]). 
For glass particle distributions with 128 neighbours, the CT 
kernel gives a density error ~ 2% larger than the CS ker- 
nel, while the HOCT4 kernel with 442 neighbours is only 
~ 1% wors«[^] (Note that the error can be very large if too 
few neighbours are used: higher order kernels require more 
neighbours to be adequately sampled.) 

In summary, formal convergence in SPH is somewhat 
subtle; it requires several important criteria to be satisfied: 

(i) increased particle number, TV; 

(ii) increased neighbour number, N n to ensure that En,j 
shrinks faster than 

(iii) a higher order kernel to maintain spatial resolution; 
and 

(iv) a kernel that is stable to clumping/banding for the 
above choice of N n . 



5 It is not clear, given these results, why |Price| j20 1(5| > argue that 
the density error is prohibitive for the CT kernel. Most likely it is 
because the error is very large when only 32 neighbours are used. 
For larger neighbour numbers, however, the performance of the 
CT and HOCT4 kernels is acceptable. 



The CT kernel with 128 neighbours and the HOCT4 kernel 
with 442 neighbours thus provide a convergent kernel pair 
that satisfy the above criteria. We will demonstrate this in 



{ 7.3 (Note that the above is simply what is required for- 
mally. It may be that for a given numerical problem Eo,» 
shrinks faster than hi without any need to raise the neigh- 
bour number. This cannot be guaranteed in general, how- 
ever.) 

The above convergence criteria - constructed simply to 
ensure that En,j shrinks faster than hi - seem rather labo- 
rious. Given that we know a priori what En,i is for each 
particle, it is tempting to simply factor it out of the mo- 
mentum equation as follows: 



dvi 
dt 



N 

£ 

3 



P^Pj 



[Pi - Pi\ ViWn 



hi pi 



-En 



(18) 



The left term is now a higher order momentum equatior|_| It 
gives zero force for constant pressure by construction, unlike 
equation [7] And, it should give much simpler convergence - 
no longer requiring increased neighbour number, or care- 



ful kernel choice (indeed we will demonstrate this in E 7.3 ) . 
However, these advantages come at a price. As pointed out 
in RHA10, notice that the left term in equation [18] is sym- 
metric in i and j inside the sum. This means that momen- 
tum is no longer conserved between particle pairs. This lack 
of manifest momentum conservation becomes a problem in 
strong shocks (Morris 1996a). We will discuss subtracted- 
En,i momentum equations and their potential for creating 
higher order SPH-like methods in a separate paper (Hayfield 
& Read 2011, in prep.). 



4 CONVERGENT FLOW & DISSIPATION 

The second problem with SPH is dealing with flow conver- 
gence - a problem common to any Lagrangian scheme. SPH 
can be thought of as both a Monte-Carlo method and a 
method of characteristics. It is a Monte-Carlo method be- 
cause a finite number of discrete particles are used to ini- 
tially sample the fluid. However, from this moment onwards 
it a method of characteristics: the particles move along 
streamlines in the flow. The problem is that the particles 
represent large unresolved patches of the fluid. Unlike real 
infinitesimal points in a fluid flow, SPH particles can ap- 
proach one another, as shown in Figure [2] This leads to 
multivalued fluid quantities at the crossing point: multival- 
ued momentum, entropy, mass, and any other fluid quantity 
that is advected with the particles. The only quantity that 
is not multivalued is the density since this is calculated by 
smoothing over a particles' nearest neighbours (see equation 



!- 



6 The left term in equation |18| is remarkably similar to the 
momentum equation discussed and proposed recently by |Abel| 
| |2011[ l. Such an equation has been proposed several times in the 
literature before (e.g. [Morris 1995 ). As was pointed out in RHA10, 
it improves mixing in SPH because it manifestly removes the En » 
error. However, if we only subtract Er^; and do nothing else, the 
method will fail because of a lack of momentum conservation in 
strong shocks, and because nothing has been done to mitigate the 
LMI. 



SPH with a higher order dissipation switch 5 




Figure 2. A schematic representation of two SPH particles ap- 
proaching one another in a convergent flow. The two particles 
carry discretely different advected quantities with them: entropy, 
Ai^2, mass rrai^, velocity Vi 2 etc. Apart from their densities that 
are manifestly smooth (c.f. equation [Cijl , all other fluid quantities 
become multi-valued at point P. Thus, we must detect when this 
situation is going to occur and add dissipation terms in all fluid 
quantities to ensure that they remain single-valued throughout 
the flow. 



The problem of particles approaching one another was 
realised very early on in the development of SPH, and led to 
the introduction of artificial viscosity. This acts to make the 
momentum between particles single valued as they approach 
one another, while maintaining energy and momentum con- 
servation. However, less appreciated in the literature is the 
need for similar dissipation terms in all other advected fluid 



quantities. This was recently highlighted by Price (20081 



For example, if two particles approach one another with very 
different entropy, their pressures will become multivalued. 
This leads to a spurious repulsive force between the parti- 
cles which inhibits mixing. In RHA10, we referred to this 
as the Local Mixing Instability (LMI). However, it can be 
thought of as a more general problem of multivalued fluid 
quantities arising when particles approach one another. 

There are two possible solutions to deal with multival- 
ued pressures in SPH. In RHA10, we used an idea from 



Ritchie & Thomas (2001) to manifestly smooth the pres- 



sures by using a modified 'RT' density estimator: 



Pi 



N 

E 



which gives (via equation! 



Pi = 



JL 1 _ 
X - 1 . "'< u '« 



(19) 



(20) 



The above density estimator, combined with equations[6]and 
[7] defines the OSPH method derived in RHA10. 

The 'RT' density estimator has the nice feature that it 
avoids multivalued pressures by construction. However, there 
is an associated cost. Consider the situation of large en- 
tropy contrasts on the kernel scale. Particles with Ai 2> Aj 
will contribute essentially zero weight, reducing the effec- 
tive kernel sampling. To maintain a constant |Eo| error, 
we must then scale the neighbour number proportional to 
the entropy contrast on the kernel scale. This becomes pro- 
hibitively expensive for astrophysically important applica- 



tions like strong blast waves. Here, OSPH gives significantly 
poorer performance than SPH for the same numerical cost. 
We demonstrate this in Appendix [B] 

For the above reasons, in this paper we take an ap- 



proach more similar to Price (20081, but with a key differ- 
ence. Price ( 2008 1 presented dissipation switches designed to 



detect (and correct) multivalued pressures. However, once 
pressures are multivalued it is already too late. As demon- 



strated recently by Valcke et al. (20101, once pressure blips 



form at flow boundaries, they cause pressure waves to prop- 
agate throughout the fluid. These damp the growth of sur- 
face instabilities and cause errors to propagate throughout 
the flow. To avoid this problem, we must detect when parti- 
cles will approach one another in advance. We can then act 
to ensure that all fluid quantities (not just the pressure) will 
be single valued by the time the particles reach one another. 
This is the strategy we adopt here. 

To detect when particles will cross, we require an accu- 
rate flow convergence detector. We take an approach similar 
tolCullen & Dehnenl (120101). iCullen & Dehnenl (120101) came 



up with the novel idea of using the time derivative of the 
velocity divergence to detect flow convergence in advance. 
They then switch on artificial viscosity to prevent particle 
inter-penetration. We use a similar idea, but consider in- 
stead the spatial derivative of the velocity divergence. As we 
will show, this has the advantage that we obtain an excellent 
estimate of the flow divergence and curl for free. 



Cullen & Dehnen ( 2010 1 focus only on the artificial vis- 



cosity. Here, we use the same switch not just for the artificial 
viscosity, but for all artificial dissipation terms. (Recall that 
we require one of these for each advected fluid quantity.) 

We have some freedom in how to construct the flow con- 
vergence detector and the dissipation terms. However, both 
must satisfy a number of constraints in order for the scheme 
to produce converged results with increasing resolution: 

(i) the switch must detect flow convergence before it oc- 
curs; 

(ii) the switch must be sufficiently robust (i.e. high order) 
as to not trigger randomly due to particle noise; 

(iii) the dissipation terms must respect conservation laws; 

(iv) the dissipation terms must shrink on a given physical 
scale with increasing resolution; and 

(v) the dissipation terms must not generate spurious pres- 
sure waves that propagate through the fluid. 

These criteria guide our choices for the switch and the arti- 
ficial dissipation terms that we describe in the following two 
subsections. The last point, in particular, is important. It is 
no good if our dissipation terms introduce more problems 
than they solve. They should act to make fluid quantities 
single valued wherever particles approach one another. But 
they should do this in a manner that respects conservation 
laws, is convergent, and does not lead to problems elsewhere 
in the flow. 



4.1 A higher order convergence detector 

We first describe our higher order flow convergence detector. 
Local flow convergence occurs wherever the velocity diver- 
gence is negative. This suggests that we should switch on 
dissipation terms if V • v, < for a given particle. However, 
if we set the magnitude of the dissipation also using V • Vj, 
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then the dissipation will only switch on once the flow is con- 
verging, not before. To detect flow convergence in advance, 
we use instead the spatial derivative of V • v; for the mag- 
nitude of our dissipation parameter ai oc ,;. This leads to the 
following dimensionless dissipation switch: 



h||V(V'Vi)| 



hflV^-Vi^ + hilV-Vil + nsCs 



V ■ Vi < 
otherwise 



(21) 



where evi OC; i describes the amount of dissipation for a given 
particle in the range [0, a max = 1]; and n s — 0.05 is a 'noise' 
parameter that determines the magnitude of velocity fluc- 
tuations that trigger the switch. Equation [21] turns on dissi- 
pation if V ■ Vi < (convergent flow) and if the magnitude 
of the spatial derivative of V • Vj is large as compared to the 
local divergence (i.e. if the flow is going to converge). 



In principle, the maximum dissipation parameter a max 
can be different for each fluid quantity. Our default in this 
paper is to use a max = 1 for all fluid variables. We investi- 
gate the sensitivity of SPHS to a ma x in Appendix [F] 



As in Cullen & Dehnen ( 2010 1, we set the local dissipa- 



tion to the above value instantaneously if cti < cn OCji : 

on — Q?ioc,i cti < aioc,; (22) 
otherwise, cti smoothly decays back to zero: 



Cti = ("loci — Qi)M Q min < Qloc.i < Cti 
Cti — ("mill — Cti)/Ti a m i n > Ct\ oc ,i 



(23) 



where Ti = /!.i/i> max .i is the timescale for the decay; and 
fm a x,i is the maximum signal velocity (Springcl 2005): 



max [i>si g ,ij] 
j 



with 



(24) 



(25) 



where Wij = — p — r 

\ r ij\ 



and Ci is the local sound speed at par- 

■ 

ticle i. 

The parameter a m i n = 0.2 ensures that the dissipation 
parameter decays all the way back to zero once particles are 
no longer converging. 

4.2 A higher order gradient estimator 



Our dissipation switch (equation 21 1 requires a good esti- 
mate of both the first and second derivatives of the veloc- 
ity field. A noisy estimator will cause the limiter to trig- 
ger unnecessarily, leading to an overly diffusive methooQ To 
achieve good quality gradients, we fit a second order poly- 
nomial to each of the fluid variables as in lMaron fc Howesl 



( 2003 1 . The first and second derivatives then follow from 
the coefficients of the polynomial fit. The full 3D algorithm 
is given in Appendix [C] Here we present a ID version to 
illustrate the idea. 



We assume that a fluid variable, qi, can be locally rep- 
resented by a smooth second order polynomial: 



qi = «o,i + a 1A Xij + a 2 , i x 2 j + 0(h A ) 



(26) 



where Xij = rij/hi. 

To determine the coefficients a n< i, we then consider the 
matrix equation Ma = q: 
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■rrijWi- 
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qjXij 
qjXij 



(27) 
(28) 



The matrix M and the vector q contain weighted moments 
that can be calculated in the usual way by summing over 
each particle's nearest neighbours. The vector a is then cal- 
culated by solving for the inverse of M. The particle gra- 
dients at the position of the particle (xtj = 0) then follow 
from q'i(0) = <zi,j and q"(0) = 2a 2 ,i- 

The above straightforwardly generalises to 3D and to 
vector fluid variables. For scalar variables in 3D we must 
solve a 10 x 10 matrix inverse to obtain a 10 coefficient fit 
(see Appendix [C|: 



2 2 

q%3 = ao,i + aijXij + 02,ij/ij + a-.i,iZij + a^jx^ + a^^y^ + 

2 

&6,iZij i 0/7 ,i^ijVij T" CLg^iXij Zij -j- ag ; iUij Zij -\- 



0(h 3 



(29) 



i/hj — [x^ , y.jj , Zjj } . 



where xy 

Note that Maron & Howes (2003) use these higher or- 
der gradients to actually move the fluid. This makes the 
method non-conservative, leading to problems in strong 
shocks. In SPHS, we use these gradients instead to conser- 
vatively maintain fluid smoothness. 

Our dissipation switch manifestly satisfies our criteria 
(i) and (ii) outlined above. It detects flow convergence in 
advance, and it is accurate since it is based on a second 
order accurate expansion of the velocity field. 

Note that a second order polynomial is the lowest order 
that we could fit in order to obtain a second derivative. In 
principle, we could fit a third or fourth order polynomial thus 
further increasing the accuracy of the switch. However, this 
comes at quite significant cost. At third order, the size of the 
moment matrix increases from 10 x 10 to 20 x 20 and requires 
an additional 40 sums over the particles to be calculated and 
stored. Secondly, for the higher order moments to actually 
help, the neighbour number should be increased. Otherwise 
noise in the third moments could make the higher order 
gradient estimator poorer than the second order estimate. 
For these reasons, we stick to the second order scheme in 
this paper. 



7 Indeed, |Rosswog| l |2010| recently advanced the idea of using 
higher order gradients for their dissipation switch. They used a 
first order accurate gradient of the pressure, whereas we use the 
gradient of the velocity divergence (which is a second derivative 
of the velocity field). 



4.3 The dissipation terms 

4-3.1 Artificial viscosity 

We start with the familiar artificial viscosity. Here, we use 



the form as in Monaghan ( 1997 1 and Springel ( 2005 ) 
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Figure 3. 3D Sod shock tube test results at t = 0.2 in SPHS (black) and SPH-CS96 (red). From left to right, the panels show: density; 
pressure; the magnitude of v x (the a;- velocity component along the shock); and the dissipation switch a (only relevant for the SPHS 
simulations). The blue line marks the analytic solution. Notice that the pressure blip at x ~ 0.2 is almost fully removed in SPHS. The 
top panels show results for a single particle-mass simulation; the bottom for the same run with multimass particles in SPHS on a uniform 
particle grid. 



Vdiss.i = - rrijllij VjWjj 



where: 

n y = 



o 



if vy ■ Tij < 
otherwise 



(30) 



(31) 



and Usig.ij and Wij are defined by 



where ay = \ [on + ay], 
equation [25] This must then generate entropy to ensure en- 
ergy conservation: 



It- 



N 

1 \ - 



2 pV ^ 



(32) 



In addition, we use a Balsara-like switch to limit viscosity in 



shear flows (Balsara||1989|. As in Cullen & Dehnen (20101, 



we apply this to our viscosity parameter ai, rather than 
directly to equation |32| This is mathematically identical, 
but means that on represents the true viscosity of the flow. 
Thus, we multiply on by a suppression function given by: 

IV -vL 



(33) 



| V ■ v\i + | V A v\i + O.OOQltv/fcj 

where c s ,i is the sound speed for particle i. Equation |33| is 
identical to the usual Balsara switch, except that we use the 
higher order gradients derived in Appendix [C] to derive the 
divergence and curl of the velocity field. 



Equations 30 and 32 satisfy our dissipation criteria (iii)- 
(v) outlined above. They respect energy and momentum 
conservation by construction; they act only on the kernel 
scale (and thus the viscosity will reduce at a given physical 
scale as the resolution is increased); and they introduce a 
numerical error only locally. 

4-3.2 Entropy dissipation 

For our dissipation in the entropy function variable Ai, we 
choose a form that explicitly conserves energy, similar to 



that proposed in Price (2008) 



N 

i 



7 — V T 

0.1 t V ■ - - hi i 1 



Ai 



7-1 



Ka (34) 



where py = [pi + pj]/2 is the symmetrised density; Kij — 
fij ■ ViWij is a symmetric smoothing kernel; Lij is a pressure 
limiter (of which more in a moment); and v^ ig t ■ is similar 
to the signal velocity, but defined to be positive definite: 



d + Cj — 3Wij if 3wij < (a + Cj) 
otherwise 



(35) 



This modified signal velocity is chosen to give more dissi- 
pation to approaching particle pairs than receding particle 
pairs. However, unlike the viscosity where the dissipation is 



fully suppressed for receding pairs (c.f. equation 31 1, we find 
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that receding pairs still require some small entropy dissipa- 
tion. This is because, while neighbouring particles can have 
discretely different velocities without serious repercussion 
(so long as they are not approaching one another), discretely 
different entropies inside the kernel will drive spurious pres- 
sure waves that affect the numerical solution everywhere. 

In fact, the above explains why adding some small en- 
tropy dissipation is preferable to doing nothing at all. The 
right amount of entropy dissipation will ensure smooth pres- 
sures and keep errors local. But the key is getting the 'right 
amount'. If we are not careful, our dissipation terms can ac- 
tually drive pressure waves and do more harm than good. 
To avoid this, we introduce a pressure limiter: 



r _ \Pi~Pj\ 



(36) 



Note that, unlike the dissipation prescription presented in 
|Price| ( [2008] ) , equation [34] poses no problem for simulations 
involving gravity. In hydrostatic equilibrium the entropy dis- 
sipation will vanish since the flow is non-converging and 
ay = 0. 

Equation 34 satisfies our dissipation criteria (iii)-(v) 



outlined above. It respects energy conservation by construc- 
tion; acts only on the kernel scale; and - through the pres- 
sure limiter - does not propagate errors non-locally. 



4-3.3 Mass dissipation (for multimass applications) 

Multimass SPH particles are very useful since they allow in- 
teresting regions of the flow to be simulated at significantly 



higher resolution (e.g. Monaghan & Vamas 1988 Meglicki 
|et al.|1993[ ). However, classic SPH runs into difficulties once 
particle masses are allowed to vary (see e.g. Ott & Schnet- 



tor 



2003 1 . The problems occur because, like the entropy, 



the particle masses are advected along with the particles. 
When particles approach one another, the masses become 
multivalued, driving a pressure wave at the mass interface. 
The problem is less severe than for the entropies because the 
masses are smoothed over in the equation of state (c.f. equa- 
tions [6] and [8| . Nonetheless, large density contrasts realised 
with multimass particles are problematic. 

Some approaches to multimass SPH have been proposed 



in the literature. Ott & Schnetter (20031 suggest adapting 



the density estimate to ensure smooth pressures by construc- 
tion - an approach very similar to the multiphase SPH pro- 



posed by Ritchie & Thomas (2001 1. Kitsionas & Whitworth 
( |2002[ ) suggest increasing the neighbour number at course- 
fine boundaries. This will act to smooth any pressure blips 
at the interface and is therefore also a viable solution. 

A key advantage of our approach is that we can treat 
any advected fluid quantity in the same manner as the en- 
tropy, above. This includes the particle masses, which allows 
us to consider a multimass SPH scheme that does not require 
raising the neighbour number at boundaries, or introducing 
a new density estimator. Treating the mass similarly to the 
entropy, above, we introduce a conservative pairwise mass 
dissipation: 



E mjj _ p 



(37) 
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Figure 4. Sod shock tube convergence test results. The x-axis 
gives the number of particles along the shock iVi£> (Figure [3] 
top, shows results for JVijj = 600); the y-axis shows the binned 
velocity error in x-bins of width 0.01; and the thick black line 
shows a scaling of N^' S (the best-possible scaling is N^). 



where = [m,+mj]/2 is the symmetrised mass. Note that 
this symmetrised mass appears only in equation |37[ and not 
in the other dissipation terms. This difference follows from 
the fact that equation [37] must respect mass conservation, 
while equations [34] and |32| - that describe the evolution of 
the specific entropy - must respect energy conservation (see 
Appendix [D] for further details). 

As for the artificial viscosity, we must then add correc- 
tion terms to ensure momentum and energy conservation. 
There is actually some freedom in how we choose to do this 
(see Appendix|D|. A simple approach is to ensure that each 
particle individually conserves its energy and momentum: 



d(rriiVi) 
dt 



= ThiVi + mii/i = 



dE l 1 . 

— — = m,u, + rriiUi + -rriiVi ■ v< + rmvi ■ v, 
at Z 



(38) 



(39) 



where we recall that Ui is the specific internal energy. Substi- 
tuting equation [38] into equation[39]then gives the correction 
terms for each particle: 



m dis 



Vdiss.z — 



; _ 1 7 — 1 rhdiss.i r 

2 pj m l 



rru 



A, 



(40) 



(41) 



We will use the above correction terms throughout this pa- 
per. (We derive a general class of correction terms in Ap- 
pendix [D] these may lead to even better results in some 
situations, but we leave this as an investigation for future 
work.) It is clear that equations |37| |40| and |4l] satisfy our 
criteria (iii) - (v) outlined above. For equal mass particle 
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applications rhdiss,; = by construction and equations |37| 
l40l and HTl have no effect. 

A final concern is that adding mass dissipation will af- 
fect our solution of the continuity equation. Taking the time 
derivative of equation [6] we have that: 

d \ ■> 



The right term is the familiar SPH continuity equation; the 
left term is a correction factor that accounts for mass dissi- 
pation. Thus, by using the familiar SPH density sum (equa- 
tion pjj), we automatically include the mass dissipation cor- 
rection to the continuity equation. However, we may still 
worry whether equation [42] tends towards equation [I] in the 
limit of infinite resolution. Substituting for rhj 
equation |42| we have: 



rn dil 



3 j,k 



)jk{m 3 - m k )K jk Wij ~ 



(43) 



where Qij = Qj t = rriij /p^okjV^^Lij , and the equation is 
very nearly vanishing since the sum is almost perfectly an- 
tisymmetric in the indices j, k (the antisymmetry is broken 
by Wij). In the continuum limit, the above term is exactly 
zero and so equation [42] does indeed tend towards equation 
[TJwith increasing resolution. 



5 TIMESTEPPING 

For our timestep control, we use individual particle 
timesteps ordered on a hierarchy of rungs in powers of two, 
as in 



Springcl (2005). Particles are placed on rungs using a 



Courant-like condition: 



(44) 



where C = 0.2 is the Courant factor. We use this same fixed 
Courant factor for all tests presented in this paper. 

In addition to the standard timestep criteria above, we 
introduce a constraint similar to that in ISaitoh fc Makino 
( |2009[ > to ensure that neighbouring particles do not differ in 
their timesteps by more than a factor of 4. 



6 IMPLEMENTATION 



We implemented our method in the GADGET- 2 code ( Springel 



2005). GADGET- 2 is a massively parallel Tree- SPH code orig- 



inally designed to model galaxy collisions, but now adapted 
for cosmological, hydrodynamic, magnetohydrodynamic and 
many other applications. SPHS acts as a new hydro module 
within GADGET-2, replacing the standard SPH parts of the 
code. We refer the reader to the original GADGET-2 paper for 
details of the gravity solver ( Springel 2005 1 . 

The SPHS hydro module, like GADGET-2 requires two 
loops over the particles. In the first loop, the densities are 
calculated (iterating to ensure equation 11 is satisfied). At 



the same time, we calculate the polynomial fluid gradients 
(for the dissipation switch; Q. In the second loop, the hy- 
drodynamic forces are evaluated along with the dissipation 



terms. Some speed comparisons between our current imple- 
mentation of SPHS and classic GADGET-2 SPH are given in 

m 



7 CODE TESTS 

In this section, we present a suite of code tests designed to 
challenge the SPHS method. In §7.2[ we use the Sod shock 
tube test to examine shocks in SPHS both with and with- 
out multimass particles, and to asses the rate of convergence 
in SPHS. In §7.3[ we use the 'Gresho' vortex test to ex- 
amine convergence in shear flows in SPHS, and the role of 
numerical viscosity. In §7.4| we use a strong Sedov- Taylor 
blast wave test to see how well SPHS performs in the pres- 
ence of extreme entropy contrasts. In §7.5[ we use a Kelvin- 
Helmholtz instability test with density contrast 1:8 - both 
with and without multimass particles - to examine mixing 
in SPHS. Finally, in fL6[ we use the 'blob' test - a 1:10 den- 
sity ratio gas sphere in a wind tunnel to assess how SPHS 
performs in more complex flow situations where shocks and 
mixing combine. 



7.1 Simulation labelling convention 

In the following subsections, we run a broad range of sim- 
ulations both in our new hydrodynamics code SPHS, and 
in 'classic' SPH (the version of SPH that is in the public 
release version of the GADGET-2 code, and that is described 
Springcl 2005). To avoid confusion, we use the following 



naming convention for these simulations: 
SPHX-KKNNx 

where X refers to the flavour of SPH: 'classic' SPH (SPH), 
or our new code (SPHS); KK refers to the kernel used (CS; 
cubic spline; equation |15[ ), (CT; core-triangle; equation |16[ ), 
(HCT; High Order Core Triangle; equation[l7]); NN refers to 
the neighbour number (42, 96, 128, 442); and x is reserved to 
describe special simulations: x = g means that the test was 
run using glass rather than lattice initial conditions; x = c 
means that the test was run using a higher order momentum 
equation (equation[l8]with the Eo term subtracted) ; and x = 
multi means that the test was run using multimass particles. 



7.2 Sod shock tube 

The Sod shock tube test is a ID tube on the interval 
[—0.5, 0.5] with a discontinuous change in properties at x — 
designed to generate a shock. The left state is described by 
pi — 1.0, Pi — 1.0, vi — 0, and the right state by p r — 0.125, 
P r — 0.1, v r = 0, where p, P and v are the density, pressure 
and velocity along the x axis. We use an adiabatic equa- 
tion of state with 7 = 5/3 and perform the test in 3D on the 
union of a 32 x 32 x 400 lattice on the left, with a 16 x 16 x 200 
lattice on the right, giving a ID resolution of Nm = 600 
points. For the SPHS simulation, we set an initial dissipa- 
tion parameter a — 1 over the initial pressure discontinuity 
(—0.05 < x < 0.05) since this test starts with a shock. We 
use lattice ICs for this test aligned with the shock, identi- 
cal to those presented in RHA10. In SPH, however, there 
is some freedom in how to lay down the particles. |Cullen| 
& Dehnen (20101, for example, use instead glass-like initial 
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Figure 5. Gresho vortex test results. The top panels show results for the HOCT4 kernel with 442 neighbours; the bottom panels for the 
CT kernel with 128 neighbours. The left panels show the v,/, velocities of particles after a time 4 = 1 for N = 64 X 64 X 8 particles; the 
analytic solution is marked in blue. The middle panels show the dissipation parameter a at the same time. The right panels show the 
three Eq error components in black, red and blue. Over-plotted is the mean binned error magnitude |Eq | (solid black line). 



conditions. These are noisier than the simple cubic lattice we 
use here, but have no preferred direction. While the choice 
of initial condition can affect the results, it should not af- 
fect the difference in the results between SPH and SPHS. 
We will demonstrate this using glass and lattice ICs for the 
Gresho vortex test in §7.3| 

The results at time t = 0.2 are shown in Figure [3] 
for SPH (red) and SPHS (black). The analytic solution is 
marked in blue. Notice that SPHS performs well on this test. 
In particular, the pressure blip at the shock, present in the 
SPH run (red), is almost completely gone. The SPH run, 
which used 96 neighbours, shows significantly more noise 
in the velocity distribution. This occurs due to symmetry 
breaking of the simple cubic lattice ICs and is reduced for 
glass-like initial conditions (Cullen & Dehnen 20101. It is 



also reduced by moving to higher order kernels that have 
larger neighbour number and are therefore correspondingly 
less noisy (e.g. Price 20101. This is why the noise is not 



present in the SPHS simulation. 

In addition, we perform this same test using multimass 



particles in SPHS, where we sample the domain uniformly 
with a 16 x 16 x 400 lattice - the same resolution as the low 
density phase in the single particle mass Sod test (Figure 
[3] bottom panels). This is an extremely challenging test for 
SPHS. The initial conditions have a sharp jump in three 
fluid variables: entropy, pressure and mass. Nonetheless, the 
solution is in excellent agreement with the analytic result. 

Finally, we perform a convergence study for the Sod 
shock tube test in Figure [4] We follow |Springei| ( |2010b| ) and 
define our error measure as: 



1 

Ll(v x ) = — — 



(45) 



where Nt is the number of bins in x, v X) i is the mean x- 
velocity in bin i, and v x (xi) is the expected analytic mean 
velocity in bin i. We use a x-bin width of 0.01. This is small 
enough to capture the improvements with increasing resolu- 
tion, but not so small as to over-sample the lowest resolution 
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Springel| < |2010b| > find that for shock tests in 2D, SPH 
performs worse than the optimal N^ scaling, giving some- 
thing closer to N^' 7 . Here we find that, by contrast, SPHS 
gives a near-optimal scaling, going as very nearly N7/ D , ex- 
cept at the highest resolution (compare the black and thick 
black lines in Figure[4]| . The slowing down of the convergence 
rate for the highest resolution simulation is due to the fun- 
damental convergence limit set by our neighbour number. 
We demonstrate this in more detail for the Gresho Vortex 
test, next. 



7.3 Gresho Vortex test 



We set up the Gresho Vortex test similarly to |Springei] 
( 2010b[ ) (and see |Gresho fc Chan|199()j ). The test involves an 
N x N x N/8 3D lattice of particles. A velocity and pressure 
field are applied to these to set up a stable vortex: 



P(R) = 




< R «: 0.2 
0.2 < R < 0.4 
R^ 0.4 



5+f R 2 
9+f R 2 - 
20i? + 41n(i?/0.2) 
3 + 4 In 2 



for OsCflg 0.2 



for 
for 



0.2 sC R s; 0.4 
i?> 0.4 



(46) 



(47) 



where R = \J x 2 + y 2 and we set p = 1 and 7 = 5/3. 

The above vortex should be stable over many rotations, 
but in practice will decay due to the numerical viscosity 
inherent in the scheme. As such, it is a useful test of the 
numerical viscosity generated in shear flows. Indeed, classic 
SPH performs poorly on this test converging very slowly to 
the wrong solution (Springcl 2010b). Such rotating config- 
urations are common in a wide class of astrophysical prob- 
lems; it is important for numerical schemes to perform well 
on such tests. 

The results for SPHS are given in Figure [5] We show, 
from left to right, the rotational velocity profile of the vortex 
(black points) as compared to the analytic solution (thick 
blue line); the dissipation parameter a, and the leading order 
error in the momentum equation Eq. We find significantly 
better performance in SPHS than was found by Springcl 
(2010b I for SPH. The primary reason for this - surprisingly 



- is not the lower viscosity of the method. The average vis- 
cosity is lower in SPHS- in the range 0.05 < a < 0.3 as 
compared to SPH that has constant a — 1 (see Figure [5j 
middle panels). But, the real reason for the improvement is 
the improved force accuracy. In Figure [5] the top three pan- 
els show the results for our default method (SPHS-HCT442), 
while the bottom three show results for SPHS using a lower 
neighbour number with the CT kernel (SPHS-CT128; see 
[7.1 for our simulation labelling convention). If anything, 
the viscosity is slightly lower for the SPHS-CT128 simula- 
tion, yet the results are worse, with increased noise and a 
bias in the rotational velocity for R ^ 0.2. The only differ- 
ence between these two simulations is the neighbour number, 
and the associated E error. Indeed, in SPHS-HCT442, the 
Eo is lower than for the SPHS-CT128 simulation (see Figure 
[5] right panels). 
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Figure 6. Gresho vortex convergence test results. The x-axis 
gives the number of particles along one side of the box Nm (Fig- 
ure[5]shows results for Nm = 64); the y-axis shows the binned 
velocity error in il-bins of width 0.01. The different line colours 
show different SPH methods; the naming convention (marked) 
is as described in j ]7.1| Three simulations: SPHS-CS42 ('classic' 
SPH); SPHS-HCT442 (our default SPHS method); and SPHS- 
HCT442e (our default method using a higher order momentum 
equation) are tested to higher resolution (Nm = 256). The thick 
black line shows a scaling of N^' 4 (the best-possible scaling is 
N~ 2 ) 



In Figure |6j we explore this further by presenting con- 
vergence tests for SPH and SPHS using varying neighbour 
number and kernel choice. We calculate the Ll(v$) error 
norm as in the Sod test (i 7.2 1, using a bin size of A_R = 0.01. 



The thick black line on the plot marks the ideal scaling of 
N^~p (ideal for a second order method away from contact 
discontinuities). The red lines show results for classic SPH 
with 42 neighbours (solid line), 96 neighbours (dotted line) 
and 128 neighbours using the CT kernel (dashed line). No- 
tice that SPH-CS42 converges very slowly with increasing 
resolution, with the error always larger than 10%. Increas- 
ing the neighbour number to 96 neighbours helps at low 
resolution but gives diminishing returns with increasing res- 
olution. This agrees with our results from RHA10, where we 
showed that in shear flows the Eo error improves only very 
slowly with increasing neighbour number for the CS ker- 
nel. This is because for neighbour number larger than ~ 40, 
the particles begin to clump preventing any significant im- 
provement in the kernel sampling. By contrast, switching 
to the CT kernel with 128 neighbours - that is manifestly 
stable to particle clumping - gives a significant improve- 
ment in the convergence rate (dashed line). Now the error 
drops to ~ 5% for Nm = 128. The solid blue line shows the 
result for SPHS using the CT kernel with 128 neighbours 
(SPHS-CT128). The results are only very slightly better 
than for classic SPH. This highlights that, for this test, it is 
the force accuracy that determines the rate of convergence, 
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Figure 7. Sedov- Taylor blast wave test results. Left three panels: The density profile of the gas at t = 0.05 for TV = 64 3 and N = 128 3 
particles for SPHS-HCT442, and for N = 128 3 particles for SPH-CS42 (using the timestep limiter described in (JS). The blue lines 
mark the analytic solution. For the SPHS simulations, the actual un-binned point particle densities are plotted in black; for the SPH 
simulation, they are plotted in grey. Notice the significantly larger noise for the SPH simulation. A mean binned profile, using a bin size 
of Ax = 0.001 is over-plotted in red. Right-most panel: Logarithmic density contours of the blast wave viewed from top at t = 0.05 
for SPHS-HCT442 with N = 128 3 particles. 



not the dissipation scheme. The solid black line shows the 
result for our default SPHS scheme: SPHS-HCT442. With 
442 neighbours and a correspondingly higher order kernel, 
our method now converges on percent level accuracy for this 
test. The convergence appears to be uninterrupted even at 
Nid = 256, though the rate is perhaps slowing. The dotted 
black line shows the results for the same simulation but run 
using glass initial conditions. The error is slightly improved, 
but the rate of convergence is identical. This demonstrates 
that our results are not sensitive to the initial particle dis- 
tribution. Finally, the solid magenta line shows results for 
our default SPHS scheme, but using a higher order momen- 
tum equation. For this simulation, we replaced equation [7J 
with the left term in equation [18] - i.e. explicitly subtract- 
ing away the E error (SPHS-HCT442e). In this case, we 
should expect to see a steady convergence rate, without any 
need to further increase the neighbour number. Indeed, this 
is what is seen. The results are now significantly better than 
any of the other SPH methods. At the very highest resolu- 
tion (A^id = 256), however, there may still be some slowing 
in the convergence rate. It is possible that this is simply a 
fluctuation (running the test at even higher resolution to 
test this seems extravagant). Alternatively, it may be that 
at these resolutions the viscosity does start to play a role, 
slowing the convergence rate. 

The most promising results for this particular test come 
from the SPHS-HCT442e method that uses a higher order 
momentum equation. Unfortunately, this same equation vi- 
olates pairwise momentum conservation between particles 
and causes problems in strong shocks (Sj3] and see Morris 
|1996a| . As such, we defer the investigation of such schemes 
to future work. We note here, however, that our default 
scheme - SPHS-HCT442 - still performs very well on this 
test, achieving percent level accuracy with increasing resolu- 
tion. As we have shown already, this is also true for the Sod 
test Q7.2 and see Figure Hb. This suggests that for most 



7.4 Sedov- Taylor blast wave 



astrophysics applications of interest, where the errors are in 
any case dominated by sub-grid physics prescriptions rather 
than the hydrodynamics solver, our default scheme - SPHS- 
HCT442 - should be sufficient, without any need to further 
raise the neighbour number. 



We set up a Sed ov- Taylor blast-wave test as in Springcl & 
Hernquist ( 2002 1 using a uniform lattice of 64 3 particles with 
We inject an explosion energy E = 1 
This corresponds to an ini- 



initial density p — 1 
into a central region r < 
tial entropy per central particle of A = 350. The remaining 
particles are assigned A — 0.05, giving an entropy contrast 
of ~ 7000. The analytic similarity solution to this problem 
is well known (see e.g. |Landau &: Lifshitz|1966[ ), and gives a 
time evolution for the blast wave radius of: 

r(t) = 1.15 ( — J (48) 

for an adiabatic index of 7 = 5/3. 

The Sedov- Taylor test is particularly challenging for 
any hydrodynamical code because of the extreme entropy 
gradient in the initial conditions. The results for SPHS for 
N = 64 3 and N = 128 3 particles are given in the left two 
panels of Figure[7] As the resolution is increased, the results 
converge on the analytic solution shown in blue: the peak 
density of the shock increases, while the low density tail 
better matches the analytic expectations. The blast wave is 
perfectly symmetric, as shown in the right-most panel. For 
comparison, the results for classic SPH (with N = 128 3 ) are 
shown in the third panel. Notice that the result is signif- 
icantly more noisy (compare the grey dots with the black 
dots in the left two panels). The reduced noise in SPHS is 
partly due to the increased neighbour number, and partly 
due to the entropy dissipation (see for example similarly less 



noisy results for this test reported in Rosswog & Price 2007 1 



The mean solution for SPH is, however, in good agreement 
with the analytic solution (compare the red and blue lines 
for the SPH-CS42 panel). Note that, for this test we had to 
use the timestep limiter described in S|5] (and see |Saitoh fc| 
Makino 2009), and so this simulation is not strictly speak- 
ing 'classic' SPH as we have defined it in this paper. Simi- 
lar results can be obtained with classic SPH by using fixed 
timesteps and a sufficiently small Courant factor (equation 
441. However, this is computationally very expensive. 

One interesting aspect of the Sedov test is that it al- 
lows us to compare the spatial resolution in SPHS-HCT442 
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with classic SPH using 42 neighbours (SPH-CS42). Notice 
that the SPH-CS42 simulation resolves higher density in 
the shock. The unbinned particles reach densities up to 
Pmax ~ 4.5 in simulation units, whereas our default scheme 
(SPH-HCT442) manages only p max = 2.7. We argued in ^3] 
that the HOCT4 kernel with 442 neighbours should degrade 
the spatial resolution by a factor / ~ 1.5 as compared to the 
CS kernel with 42 neighbours. Since the shock front for the 
Sedov test is one dimensional, then we can expect a lower 
peak density in SPHS-HCT442 of a factor ~ /. This is al- 
most exactly what is seen since p m ax,SPH/p m ax,SPHS = 1-67. 
In practice, however, the spatial resolution of the SPH sim- 
ulation is not this good because of the increased noise. We 
ought to trust only the averaged solution that is shown in 
red. For this averaged density, the peak is significantly lower, 
with pmax ~ 3.3 - only a factor 1.2 better than our default 
SPHS method. We conclude from this that SPHS-HCT442 
does not significantly degrade the spatial resolution as com- 
pared to SPH-CS42, especially once the reduced noise of the 
method is taken into account. 



7.5 Kelvin-Helmholtz test 

We set up a 1:8 density contrast Kelvin-Helmholtz (KH) test 
in 3D as in RHA10. We used a periodic thin slab defined by 
x e {-0.5, 0.5}, y 6 {-0.5, 0.5} and z G {-1/64, 1/64}. The 
domain satisfied: 



P, T, v x = 



pi,T\,v\ \y\ < 0.25 
P2,T 2 ,v 2 M>0.25 



(49) 



The density and temperature ratio were R p = p\j pi = 
T2/T1 = c|/ci, ensuring that the whole system was in pres- 
sure equilibrium. The two layers were given constant and op- 
posing shearing velocities, with the low density layer moving 
at a Mach number M2 = —V2/C2 ~ 0.11 and the dense layer 
moving at Aii = Ai 2 \/Rp- The density ratios considered in 
this work are small which assures a subsonic regime where 
the growth of instabilities can be treated using equation |52| 
( Vietri et al.|1997| >. 

To trigger instabilities, velocity perturbations were im- 
posed on the two boundaries of the form: 



5v y [sm(2n(x + A/2)/A) exp(-(10(y - 0.25V/) 
- sin(27nr/A) exp(-(10(y + 0.25)) 2 )] 



(50) 

where the perturbation velocity Sv y /v = 1/8 and A = 0.5 is 
the wavelength of the mode. 



The linear growth rate of the KHI is given by: (Chan 



drasekhar 1961 ) 



w — k 



(P1P2) 1/2 V 
(Pi + P2) ' 



(51) 



where k = 27r/A is the wavenumber of the instability, pi and 
p 2 are the densities of the respective layers and v — v\ — V2 
is the relative shear velocity. The characteristic growth time 
for the KHI is then: 



TKH 



2_7T _ (pi + p2)\ 
W {plp2) 1/2 v' 



(52) 



We set up two simulations to satisfy the setup de- 
scribed above. An equal mass particle simulation with N = 



2, 359, 296, and a multimass version with N = 524, 288. The 
latter simulation used a uniform grid of particles, with mass 
ratio 1:8 to describe the density step. To satisfy pressure 
equilibrium everywhere, the temperatures were adjusted at 
the boundary to be consistent with the SPH density step 
that is smooth (c.f. equation |6J. 

The results of the test at times tkh = 1,2 and 3 
are shown in Figure [8] The top row shows the results for 
classic SPH (SPH-CS42); the middle row shows the results 
for our default SPHS scheme (SPHS-HCT442), using sin- 
gle mass particles; the bottom row shows the results for 
SPHS-HCT442 using multimass particles. The SPH results 
are poor, with no mixing observed between the fluid layers, 
similarly to what has been reported in previous works (e.g. 



Agertz et al.||2007| RHA10). By contrast, the SPHS results 
show the growth of KH rolls on the correct timescale and re- 
solved mixing into the fully non-linear regime. Furthermore, 
the single mass simulation and multimass simulation (mid- 
dle and bottom panels) are in excellent agreement. They dif- 
fer in the details of the non-linear evolution caused by the 
growth of smaller noise-seeded rolls. But the growth time for 
the primary KH roll is in excellent agreement with analytic 
expectations, while the non-linear evolution is qualitatively 
similar. The multimass simulation is slightly more diffusive 
due to the additional mass dissipation between particles at 
the boundary. However, this simulation (because of the lower 
particle number) ran almost 5 times faster. 

As discussed in our previous paper (RHA10), the im- 
proved performance for the KH test in SPHS is a result 
of both the improved force accuracy (due to the increased 
neighbour number and higher order, stable, kernel), and the 
improved dissipation. We demonstrate this in Appendix |E| 
where we show the effect of switching off the entropy and 



mass dissipation terms (equations 34 and 37 1 for this test 



7.6 The 'blob' test 

The 'blob' test is a spherical cloud of gas of radius R c \ in 
a wind tunnel with periodic boundary conditions. The am- 
bient medium is ten times hotter and ten times less dense 
than the cloud so that the system is in pressure equilib- 
rium. The wind velocity (u w ind = c s M) has an associated 
Mach number M = 2.7. This leads to the formation of a 
bow shock after which the post-shock subsonic flow inter- 
acts with the cloud and turns supersonic as it flows past it. 
The test was first presented (with a full analytic analysis) in 
~| 2007) , 



Agertz et al. 



Here, we set up the test as in RHA10 
with N = 126, 744 in the blob, arranged on a lattice. As in 
RHA10, we seed an initial inwards perturbation on the blob 
surface. 

The results at times tkh = 1 (top), 2 (middle), and 3 
(bottom) in classic SPH (left) and SPHS (right) are given 
in Figure [9] In classic SPH, similarly to what has been re- 
ported in previous works, the blob does not break up and 
survives for the full length of the simulation. Furthermore, 
the suppression of instabilities at the fluid interface is suf- 
ficient to remove the inward perturbation that was seeded 
in the initial conditions. By contrast, this perturbation is 
clearly visible in the SPHS simulation and grows causing 
the blob to split down the middle in good agreement with 
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Figure 8. KH1:8 test results. From left to right panels show: logarithmic density contours at: tkh = 1,2 and 3; and the dissipation 
parameter a at tkh = 1- The top panels show results for a single particle mass simulation (TV = 2,359,296); the bottom panels for a 
multimass simulation with a uniform density particle distribution (TV = 524,288). 



both Eulcrian codes and our OSPH method (RHA10). Fi- 
nally, notice that the symmetry of the blob is well-preserved 
even at tkh — 3 - well into the non-linear regime. 



tial conditions (to avoid a divide by zero in the percentage 
errors). The worst performance is for the Sedov- Taylor test 
that conserves energy at the ~ 5% level. However all other 
tests conserve energy, mass, momentum and angular mo- 
mentum to better than 1% over the full simulation time. 



7.7 Conservation 

Figure [10] summarises the conservation performance of 
SPHS for all of our tests. From left to right, we show the 
conservation of energy, momentum, angular momentum and 
(where relevant) mass. The results are normalised to a sim- 
ulation time of 1, where "1" is the maximum time presented 
in this paper (i.e. for the KH1:8 test this is tkh = 3). Mo- 
mentum and angular momentum conservation results are 
only shown where these quantities are not zero in the ini- 



7.8 Code performance 

Figure [TT] compares the ratio of the speed of our default 
SPHS scheme (SPHS-HCT442; black squares) and SPH- 
CS442 (red squares) to classic SPH (SPH-CS42). (We use 
the simulation naming convention as described in |7.1| ) For 
all tests, we used 16 processors. The Sod tests were compared 
at the -/Vid = 600 resolution; the vortex tests at TVid = 128; 
and the Sedov tests at TVid = 128. In all cases, we eliminate 
the start-up time costs (the time taken to complete step 
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Figure 9. Blob test results at tkh = 1 (top), 2 (middle), and 3 (bottom) for classic SPH (left) and SPHS (right). All plots show 
logarithmic density contours. 




Figure 10. Conservation in SPHS. From left to right, the panels show conservation of energy, momentum, angular momentum, and 
(where relevant) mass for the simulation test suite presented in this paper. The coloured lines show results for: the multimass Sod test 
(black); the Sedov- Taylor test (red; for TV = 128 3 particles); the Gresho Vortex test (green; for N = 64 X 64 X 8); the multimass KH1:8 
test (blue); and the Blob test (purple). The results are normalised to a simulation time of 1, where "1" is the maximum time presented 
in this paper (i.e. for the KH1:8 test this is tkh = 3). Momentum and angular momentum conservation results are only shown where 
these quantities are not zero in the initial conditions (to avoid a divide by zero in the percentage errors). 



zero). There is some significant variation in speed across all 
of the tests with the cost of SPHS ranging from 2 to 4 times 
that of SPH-CS42, but typically SPHS is 3-4 times slower 
at like particle number. 

Note that the above speed tests are conservative. We 
could equally well conduct the tests at like numerical er- 



ror, rather than like particle number. For the Gresho vor- 
tex test, for example, it is unlikely that SPH-CS42 can ever 
achieve equivalent accuracy to SPHS for any reasonable par- 
ticle number (see Figure|6|. To obtain ~ 1% accuracy on this 
test, classic SPH would require an enormous particle number 
and be significantly slower than SPHS. 
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Figure 11. The ratio of the speed of our default SPHS scheme 
(SPHS-HCT442; black squares) and SPH-CS442 (red squares) to 
'classic' SPH (SPH-CS42). (We use the simulation naming con- 
vention as described in |7.1| ) For all tests, we used 16 processors. 
The Sod tests were compared at the N±j} = 600 resolution; the 
vortex tests at Nny = 128; and the Sedov tests at iVirj = 128. In 
all cases, we eliminate the start-up time costs (the time taken to 
complete step zero). 



Finally, we have not made any attempt to optimise our 
current implementation of SPHS. Faster neighbour search 
algorithms, or neighbour caching could conceivably gain 
back much of the speed losses as compared to classic SPH. 
In addition, for real astrophysics applications, the additional 
work done on the neighbours may be compensated by im- 
proved timestepping (due to the reduced noise), and better 
load balancing in highly parallel simulations. Such consider- 
ations are beyond the scope of this present work. 



8 CONCLUSIONS 

We have presented an implementation of Smoothed Particle 
Hydrodynamics (SPHS) that has two novel features. The 
first is an improved treatment of dissipation. We use the 
spatial derivative of the velocity divergence as a higher or- 
der dissipation switch. Our switch - which is second order 
accurate - detects flow convergence before it occurs. If par- 
ticle trajectories are going to cross, we switch on the usual 
SPH artificial viscosity, as well as conservative dissipation 
in all advected fluid quantities (for example, the entropy). 
The viscosity and dissipation terms (that are numerical er- 
rors) are designed to ensure that all fluid quantities remain 
single-valued as particles approach one another, to respect 
conservation laws, and to vanish on a given physical scale 
as the resolution is increased. The second novel feature is 
the use of significantly larger neighbour number (442) to 



improve the force accuracy. As in our previous work, we use 
a novel kernel function that is: (i) higher order such that 
the spatial resolution is not significantly degraded by our 
larger neighbour number; and (ii) that has a constant first 
derivative in the centre to prevent particle clumping (this 
latter ensures a smooth particle distribution on the kernel 
scale, which is necessary to obtain the improvement in the 
force accuracy). 

We have demonstrated that SPHS alleviates a number 
of known problems with 'classic' SPtQ successfully resolv- 
ing mixing, and recovering numerical convergence with in- 
creasing resolution. An additional key advantage is that - 
treating the particle mass similarly to the entropy - we are 
able to use multimass particles, giving significantly improved 
control over the refinement strategy. 

We have presented a wide range of code tests: the Sod 
shock tube, Sedov- Taylor blast wave, Gresho vortex, Kelvin- 
Helmholtz instability, the 'blob test', and some convergence 
tests. Our method performed well on all tests, giving good 
agreement with analytic expectations. For some tests, like 
the Gresho vortex, most of the improvement over 'classic' 
SPH is due to the improved force accuracy. For other tests 
like the (high density contrast) Kelvin-Helmholtz instability, 
the improved dissipation is equally important. We deliber- 
ately picked challenging tests that involve sharp features in 
one or more of the fluid quantities. These are inherently dif- 
ficult to resolve for our method that is manifestly smooth, 
yet we show that SPHS copes well even in such situations. 

In our current implementation (that is likely sub- 
optimal) SPHS is typically 3-4 times slower than 'classic' 
SPH (using 42 nearest neighbours) for like particle number. 
However, this additional cost should be offset against the 
improvement in the quality of the hydrodynamic solution 
in SPHS, the significantly reduced noise, and the improved 
rate of convergence. For the Gresho vortex test, for exam- 
ple, SPHS achieves ~ percent level accuracy as compared to 
~ 10% in SPH for the same particle number. 

The main remaining flaw in the SPHS algorithm is its 
low order. This means that formal convergence requires the 
neighbour number to be increased along with the particle 
number (using increasingly higher order stable kernels to 
maintain spatial resolution). However, our default kernel 
choice with 442 neighbours is already sufficient to obtain 
percent level accuracy on the hydrodynamic tests we present 
here. It is unlikely that the neighbour number will need to 
be increased further than this for most astrophysical appli- 
cations of interest. 

SPHS will be useful for any astrophysics application 
involving multiphase fluid flow (e.g. resolving the ISM in 
galaxy discs), or where the use of multimass particles would 
be advantageous. We will apply it to a broad range of prob- 
lems in forthcoming papers. 
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APPENDIX A: A FULLY CONSERVATIVE 
VERSION OF SPHS 

A fully conservative version of SPHS can be constructed 
by replacing equation [7] with that in S pringel fc Hernquistj 
(20021): 



JV 

— — = — > m, 

3 



(Al) 



where the function fi is a correction factor that ensures en- 
ergy conservation for varying smoothing lengths: 



fi= 1 



hi dpi 
3pi dhi 



(A2) 



As discussed in RHA10, the above momentum equation 
gives improved (in fact manifest) energy conservation, but 
larger truncation error. For applications where energy con- 
servation is of paramount importance (for example, where 
a system is evolved for many dynamical times), the above 
equation should be used. However, in this case, care must 



also be taken over the timestepping (e.g. Dehnen & Read 
2011 1. For the tests presented in this paper, the energy losses 



due to variable timesteps dominate and the above momen- 
tum equation gains only ~ 0.5% in energy conservation for 
the KH1:8 multimass test (see \ 7.7 1; and ~ 2% for the Sedov 
test (a factor ~ 2 improvement in both cases). However, the 
larger truncation error introduces significantly more diffu- 
sion. This is shown in Figure [AT] where we present results 
for the multimass KH1:8 test (see j |7.5| l at time tkh = 1, 
using equation |A1| For this reason, our default choice for 
SPHS is the momentum equation [7] 



APPENDIX B: 
PRESSURES 



THE TROUBLE WITH 'RT' 



In RHA10, we used the same equations of motion as SPHS, 
but with the pressure estimator in equation |20| This ensured 
manifestly smooth pressures throughout the flow, allow- 
ing us to successfully model mixing between different fluid 
phases. However, while equation |20] gives excellent perfor- 
mance for multiphase flow applications, it performs poorly 
in strong shocks where the entropy gradients on the ker- 
nel scale are large. We show this in Figure |B1| where we 
plot results for the Sedo v- Taylor blast wave problem (with 



N = 128 ; and see [7.4|, using the 'RT' pressure estimator 
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Figure Al. KH1:8 multimass test results at tkh = 1 for a fully 
conservative version of SPHS (that uses equation |A1| instead of 
equation[7]l . Notice that the results are significantly more diffusive 
than our default scheme shown in Figure [8] 



the resulting shock front, while very smooth, is not in good 
agreement with the analytic curve shown in blue. 

The above highlights the key problem with 'RT' densi- 
ties and pressures. Particles in the kernel with very different 
entropies are down-weighted in the sum. This means that 
to obtain good kernel sampling, we must scale the neigh- 
bour number with the entropy contrast on the kernel scale. 
For the Sedov- Taylor blast wave, where the initial entropy 
contrast is ~ 7000, this is prohibitively expensive. Not do- 
ing this, however, leads to a significant numerical error as 
can be seen in Figure |BT) For this reason, in this paper, we 
have abandoned the density and pressure estimators given 
in equations 1 1 9| and |20| Instead, we ensure smooth pressures 
through our higher order dissipation switch described in Q 



APPENDIX C: FITTING AN 7VTH ORDER 
POLYNOMIAL TO A FLUID QUANTITY 

We describe here an algorithm for fitting an order N poly- 
nomial to an irregular point distribution (see e.g. |Fan fc| 
Gijbels|[l996l |Maron fc Howes||2003[). We give the relevant 



equations for a second order fit in three dimensions, but the 
method straightforwardly generalises to arbitrary order and 
dimension. Assuming that a fluid quantity qi defined at par- 
ticle position i is smooth (and therefore differentiable), we 
can perform a second order polynomial expansion at a point 
j about i: 



q%i = ao,i + ai,iXij + a 2 ,iyij + 03,iZij + a^ititj + as,iVij + 

2 



0(h 3 



(CI) 



(equation 20 1 without entropy dissipation. As can be seen, where xy 



j/hi 



[Xij , Vij i Zij\ 
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~ + is the symmetrised smooth- 



arid XV i_ 

ing kernel (the superscript T means transpose). 

Having determined all of the coefficients of a (by solving 
a = M _1 q), the gradients of q evaluated at i then simply 
follow as: 

dqt dqi dqt 



= a 2 \ 



= a 3 



dx 1 dy " dz 
and similarly for the second derivatives. 



(C6) 



Figure Bl. Sedov- Taylor blast wave test results using the 'RT' 
pressure estimator (equation |20[ | without entropy dissipation. The 
plot shows the density profile of the gas at time t = 0.05, similarly 
to Figure [7] Notice that the shock front is displaced with respect 
to the analytic curve (blue). 



The coefficients of this expansion can then be deter- 
mined by inverting the following 10 x 10 matrix equation: 



Ma = q 

where: 



[do, ax, 0,2, 0,3, 04, 05, as, 0-7, as, ag\ 



(C2) 



(C3) 



APPENDIX D: A GENERAL DERIVATION OF 
CONSERVATION TERMS FOR MULTIMASS 
SPHS 



In section pi. 3. 3| we introduced a multimass dissipation term 
for SPHS. This requires some correction terms to restore 
energy and momentum conservation. In this appendix, we 
derive the general class of such correction terms. 

Our dissipation terms must obey mass, momentum and 
energy conservation: 



M = = ^2 1 



MV = = ^ m i v J = mjirj 



+ m,v. 
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+ Uj)+ mj (vj ■ Vj + Uj) (D3) 



First, let us verify that equation [37] satisfies equation |D1 
Substituting for ihj — mdiss.j, we have: 



M = Qjk(mj - m k )K jk = 

3,k 



(D4) 



where Qi 



Qj 



rr%ij I PijOLijvl: sii Lij and the above is 



(%3 — Vji — 11H3I PijU-13 u sis,ij % 3 

zero because it is antisymmetric in j, k. Note that this ex- 
plains why we must use a symmetrised mass in equation |37| 
Qij must be symmetric in order to ensure mass conservation. 
Now, let us substitute rhj = mdi ss ,j into equation |D2| 



= ^2 m jVdiss 



J T ^ ^Cj 

J,k 



k (nij - m k )K jk Vj 



(D5) 



where we have split the acceleration into a dissipative cor- 
rection term, and all other normal SPHS terms: V; = 
Vdiss.i +v rcst ,i, and then used the fact that . mj v rcs t,j = 
by construction for SPHS. 

We may now select any form we like for Vdi ss ,i so long 
as it satisfies equation |D5| In §4.3.3[ we chose a form that 
conserves momentum on a per particle basis, but we may 
also choose a form that fluxes the momenta, for example: 
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(D6) 



It is straightforward to show that the above correction term 
also conserves momentum since it makes equation [D5] anti- 
symmetric in j, k. 

We may then derive a similar constraint equation for 
our energy correction term. As an example, let us substitute 
rhj = rhdiss.j and equation ID6I into equation ID3I This gives: 



= ^ {Qjk(rrij - m k )K jk ) 

(m-j - m k ) 
li k ~ 



3,k 



+ 



KjkVk + UdissJ 



(D7) 



where similarly to the above, we have dropped all contri- 
butions involving the standard SPHS terms since these are 
already conservative and therefore vanish. 

We may then derive a correction term for 



3 



-Ki 



l 



V,- + Uj 



(D8) 



It is straightforward to verify that substituting equation |D8| 
into equation |D7| makes the equation antisymmetric in j, k 
and thus restores energy conservation. 

It is clear from the above examples that we may use the 
above constraints to derive a whole class of correction terms. 
Some of these may give better performance than equations 
|38| and [39| that we use as our default in this paper. Such a 
study is, however, beyond the scope of this present work. 



APPENDIX E: THE IMPORTANCE OF 
DISSIPATION TERMS FOR MULTIPHASE AND 
MULTIMASS FLOW IN SPHS 

In this Appendix, we show the effect of switching off our dis- 
sipation terms in entropy and mass for the KH1:8 multimass 
test ({ 7.5 1. The results are shown in Figure \ 



El 



As expected, 

the entropy dissipation is extremely important: without it 
there is no mixing between the different fluid phases (see left 
panels). Notice, however, that even without dissipation, the 
KH rolls do grow on the correct timescale unlike in the clas- 
sic SPH simulation (Figure [8] top row). This demonstrates 
(similarly to our findings in RHA10) that the improved force 
accuracy in SPHS is responsible for the correct growth rate 
of the rolls, while the improved dissipation is responsible for 
actual mixing between the different fluid phases. 

The effect of the mass dissipation is more subtle. With- 
out mass dissipation, mixing is also inhibited, but the effects 
are less strong than for the case without entropy dissipation 
because, unlike the entropies, the masses are smoothed in- 
side the density surrj^] (equation [fij). 



9 In fact, the results for this test without mass dissipation are 
rather similar to the KH1:8 single mass test we presented using 
OSPH in RHA10. This similarity arises because in OSPH the 
entropy — like the particle masses — is smoothed inside the pressure 
estimator (equation |20|. 



APPENDIX F: THE SENSITIVITY OF SPHS 
TO THE DISSIPATION PARAMETERS 

In this appendix, we assess how sensitive SPHS is to the 
choice of dissipation parameters. As our default, we have as- 
sumed a single dissipation parameter for viscosity, mass dis- 
sipation and entropy dissipation: a = a v = a m = a a ~ 1. 
This default choice is natural from the definition of the dis- 



sipation/viscosity equations 30 32 34 and 37 These assert 
that the dissipation should proceed proportional to the jump 
in the given fluid quantity (mass, entropy etc.) and on a 
timescale set by the signal velocity. Thus, we expect a nor- 
malisation parameter in each case of order unity. Nonethe- 
less, a is a free parameter and we should check that our 
results are not sensitive to it. To test this, in Figure |FT) we 
consider the effect of varying a v and ola for the Sod shock 



tube test ({ 7.2 1 at two different resolutions. 



From Figure |F1| we see that our results are not sensi- 
tive to the entropy dissipation parameter a a- Over a wide 
range 0.1 < ola < 5, the results change only slightly. More 
importantly, the differences decrease with increasing resolu- 
tion (compare the red and black lines in the left two pan- 
els of Figure [FlJ. The results are more sensitive, however, 
to the choice of viscosity parameter a v . For low viscosity 
{piv = 0.1), we have spurious oscillations in the solution. 
Reassuringly, however, these decrease with increasing reso- 
lution (compare the green and black curves in the right two 
panels of Figure Fl I. The results are poor, however, if a v is 
too large. For a v = 5, there is a strong under-shoot in the 
density at the shock that does not improve with increasing 
resolution. 

We conclude that the results in SPHS converge with 
increasing resolution independently of a v or a a, so long as 
a v is not too large. 
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Figure El. KH1:8 multimass test results: the effect of removing the dissipation terms. From left to right the panels show: logarithmic 
density contours at: tkh = 1 (top) and tkh = 2 (bottom) for SPHS run without entropy or mass dissipation (left); without mass 
dissipation (middle); and the full SPHS scheme (right). (The right panels reproduce the results from Figure [s]) 
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Figure Fl. Sensitivity to the dissipation parameters for the Sod shock tube test (j ]7.2| ). From left to right, the panels show the density 
profile at t = 0.2 (similarly to Figure p5|, for varying entropy function dissipation parameter a^, and viscous dissipation parameter a v , 
at low resolution Nip = 200 and higher resolution Nijj = 400, as marked. The blue line marks the analytic solution. Notice that so long 
as a v is not too large, the results converge with increasing resolution independently of the choice of dissipation parameters. 
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